A study of a diauxic growth experiment using an expanded dynamic flux balance framework

Flux balance analysis (FBA) remains one of the most used methods for modeling the entirety of cellular metabolism, and a range of applications and extensions based on the FBA framework have been generated. Dynamic flux balance analysis (dFBA), the expansion of FBA into the time domain, still has issues regarding accessibility limiting its widespread adoption and application, such as a lack of a consistently rigid formalism and tools that can be applied without expert knowledge. Recent work has combined dFBA with enzyme-constrained flux balance analysis (decFBA), which has been shown to greatly improve accuracy in the comparison of computational simulations and experimental data, but such approaches generally do not take into account the fact that altering the enzyme composition of a cell is not an instantaneous process. Here, we have developed a decFBA method that explicitly takes enzyme change constraints (ecc) into account, decFBAecc. The resulting software is a simple yet flexible framework for using genome-scale metabolic modeling for simulations in the time domain that has full interoperability with the COBRA Toolbox 3.0. To assess the quality of the computational predictions of decFBAecc, we conducted a diauxic growth fermentation experiment with Escherichia coli BW25113 in glucose minimal M9 medium. The comparison of experimental data with dFBA, decFBA and decFBAecc predictions demonstrates how systematic analyses within a fixed constraint-based framework can aid the study of model parameters. Finally, in explaining experimentally observed phenotypes, our computational analysis demonstrates the importance of non-linear dependence of exchange fluxes on medium metabolite concentrations and the non-instantaneous change in enzyme composition, effects of which have not previously been accounted for in constraint-based analysis.


Introduction
Computer models are invaluable tools in capturing and systematizing new knowledge, especially for the complex phenomena found in biology. Genome-scale metabolic models (GEMs) contain parts of a solution [16,20], thus requiring decisions to be made on the part of the modeler that may produce different results from the same underlying data. Such decisions are always a part of modeling, but in the case of dFBA, they are not usually explicit, leading to a number of different solutions to the same basic problems. A notable exception is the modeling framework COMETS, which hosts a range of functionalities and provides interfaces for several popular programming languages to its open-source code, although most of that code is implemented in Java [18,21]. Overall, the tools for performing dFBA are generally aimed at expert users who are either able to improvise the required coding themselves or have very specific aims.
Most biochemical reactions need to be catalyzed by enzymes to proceed at a physiologically relevant rate (or at all). As an enzyme both has a finite rate of catalysis and a mass, this places a constraint on the possible productivity of a given amount of cells. While there are a few different approaches to incorporate this as a limitation in FBA [9,13,22], they all spring from the above assertions and assumption that this carries additional, meaningful constraints on the flux distribution in a metabolic network. Here, we collectively refer to these approaches as enzyme-constrained FBA (ecFBA). Besides the sound arguments underlying the ecFBA approaches, they also improve the ability of the constraint-based framework to reproduce some phenomena observed experimentally that are outside the capability of basic FBA without the addition of extensive and somewhat arbitrary constraints. An example is the modeling of overflow metabolism, where glucose is incompletely metabolized, thus yielding a much lower amount of energy per glucose molecule, yet a higher amount of energy per time [9,23]. While the ecFBA approaches are useful in many situations, these models are data-hungry and require extensive information on enzymes associated with different reactions, namely their mass and the turnover numbers [9]. As these enzyme-constrained models can be somewhat cumbersome to build, requiring extensive retrieval processes from databases such as SABIO-RK [24] and BRENDA [25], efforts have been made to facilitate the process through automation and more standardized protocols [9,14,26].
A common application of FBA is to study the interaction between a given organism and its environment over time, such as in the case of a batch fermentation [11]. While it is not strictly necessary to run a dFBA simulation of the entire process to do so, batch fermentation provides a natural testing ground for assessing the accuracy of both the GEM and the model of the fermentor itself, as well as testing the validity of the different assumptions at work. It is also a display of a GEMs predictive powers that might be easier to grasp than more domain-specific analysis methods such as phenotype-phase-plane analysis [27]. Due to the potential profits of improving the efficiency of industrial batch production of added-value bio-molecules, and the requirements for model fidelity in order to do so, a natural step then becomes to combine the simpler-to-interpret aspects of dFBA with the increased fidelity of ecFBA [28][29][30][31], a combination we will call decFBA in the following. The DOA approach requires added theoretical as it increases the complexity of both formulating and solving the mathematical problems involved. In the case of SOA, it can be implemented in a straightforward way. One of the key merits of ecFBA is that it improves a GEM's ability to reproduce certain metabolic phenomena crucial to the accurate description of cellular behavior, such as overflow metabolism [9].
Also required for accurate simulation is the BOF, which is defined as a pseudo-reaction, listing the amount of different (pseudo-) metabolites needed per hour to produce a growth rate of 1 h −1 . As a model of this type operates in a gCDW −1 capacity, this becomes the transient amount of different metabolites needed for 1 gCDW to produce 1 gCDW h −1 , without explicitly considering the compounding effects of exponential growth [1]. As the BOF serves a central role in GEMs, determining the overall flux distribution by being the objective of the entire system it is trying to optimize for, it should receive special care and attention. While the traditional approach has been to approximate it based on data compiled from literature, often from different strains or conditions, the subject has been the target of increasing interest in recent years [2][3][4][5][6]. Yet, much work remains to be done, both in terms of establishing protocols and analysis techniques [6], and in terms of a theoretical framework for integrating the resulting data [5] before we can reasonably expect BOFs to allow accurate predictions in varying and novel environments. Due to its central role in metabolism, the BOF is of special importance when trying to quantitatively predict metabolic behavior. Thus, the simulation of a concrete situation can be a useful test for verifying a given GEM biomass composition. This is especially true when looking at parameters such as terminal biomass, where the amount of biomass produced in a given fermentation should be predictable in a simulation using a GEM containing an accurate BOF.
In the following, we expand decFBA using the SOA approach by explicitly adding constraints on how rapidly the mass fraction occupied by enzymes can be reallocated, named enzyme change constraints (ecc). We demonstrate how this has significant effects on modeling dynamics by showing its impact on GEM behavior during simulation of generic diauxic growth on glucose. Finally, we proceed to test the impact of the ecc on a GEM's ability to recreate a specific fermentation experiment, and how the ecc can help to explain observed behavior.

Batch fermentation
A carbon-limited batch fermentation of E. coli BW25113 (from the Keio collection [32]) (hereafter E. coli), was carried out in a 3 L Eppendorf NewBrunswik BioFlo 115 bioreactor in batch setup, using minimal glucose media as specified in S1 File (the NaCl was supplied by VWR, while the rest of the compounds were supplied by Sigma-Aldrich). Glycerol stock solution samples prepared from a single colony were grown in 250 mL baffled shake flasks overnight at 37˚C and 200 rpm shaking in glucose minimal M9 medium. The fermentors containing 1.5 L of sterile medium were then inoculated with 40 mL of culture. The fermentation parameters were 37˚C, pH 7, and 40% dissolved oxygen (DO). During setup, the pH electrodes were calibrated with pre-mixed solutions of pH 4 and 7. The pH was prevented from dropping below 7 through automated addition of NaOH (4 M). The DO electrodes were calibrated to 0% by flushing the electrode for 20 min with nitrogen gas, and to 100% by running the reactor with the medium and filtered air inflow of � 600 mL min −1 at 500 rpm stirring until equilibrium. Stirring was then coupled to DO, ensuring an oxygen saturation in the medium of � 40%, with 200 rpm and 1000 rpm as lower and upper bounds. Analysis of the off-gas composition and flow-rate was performed using a ThermoScientific Prima BT Benchtop MS, which was calibrated before each run.
Sampling was performed every 30 minutes for optical density (OD) measurements and medium analysis, and every 30-60 minutes for dry weight. The OD was measured using a spectrophotometer at a wavelength of 600 nm by sampling 5 mL. An aliquot of the sample was diluted with purified water such that the measured OD was between 0.2 and 0.6. The medium samples were taken as 2 mL aliquots from the OD samples. They were spun at 6700 RCF for 2 min, the supernatant was sterile filtered (0.2 μm polyethylensulfon syringe filters) and stored at −20˚C until it was analyzed. For dry weight measurements 10 mL were sampled, washed (centrifuged and redissolved) once with 5 mL 0.9% NaCl before being centrifuged again and redissolved in purified ion-free water, and added to a dried and weighed aluminum pan. The pans were completely dried (for � 36 h at 120˚C) before being weighed again.

Medium analysis
The quantification of medium constituent concentrations was performed by NMR, using ERE-TIC2 [33] in the Bruker TopSpin 4.0.8 software. The protocol is based on Søgaard et al. [34] and was further developed.
The samples were thawed and mixed 1:10 D 2 O with 0.75% TSP (purchased from Sigma), before 600 μL of the resulting solution was added to a 5 mm 7 inch NMR tube. It was then analyzed using a 600 MHz (14.1 T) Bruker NMR spectrometer running with the H 2 O + 10% D 2 O solvent setting and water suppression ( 1 H NMR, noesygppr1d). The acquisition parameters were: 4 dummy scans, 32 scans, SW 20.8287 ppm, O1 2820.61 Hz, TD 65536, TE 298.0 K, D1 4 s, AQ 2.621 439 9 s. P1 was calibrated for each sample to ensure accurate quantification. A creatine solution (70 mM in D 2 O-TSP solution) was used as an external standard, quantified using the singlet at � 3 ppm. Compounds responsible for peaks in the spectra were identified based on reference 1 H-NMR spectra available in the Human Metabolome Database (HMDB) [35][36][37], as well as the software Chenomx [38] and published literature by Fan [39]. Note that, glucose quantification was not based on the α-Glucose doublet at � 5.2 ppm, as these peaks can be affected by the water suppression sequence [40]. Instead, we used the quartet at � 3.2 ppm for this quantification (reference HMDB compound ID HMDB0000122).

GEM, COBRA and code availability
The GEM used for all modeling was the E. coli iJO1366 [41] model, or more specifically, its enzyme-constrained modification iJO1366 � as presented by Bekiaris [14]. The coding and modeling was performed with MATLAB 2020a, using version 3.0 of the COBRA Toolbox for MATLAB [20], and the linear program solver for version 2.7 of the GNU Linear Programming Kit (glpk) solver (https://www.gnu.org/software/glpk/). Using the glpk solver, the feasTol parameter had to be lowered from the default 1e − 6 to 1e − 9 to prevent numerical issues with the algorithm placing limits on enzyme reallocation.
The code produced and used in the paper is available in S3 File, and on GitHub, at https:// github.com/emikar/FermentationSim. Anyone is permitted and encouraged to freely copy the code and implement their own modifications as needed. Details for the software can be found in S2 File.

Simulating the fermentation
The batch fermentation performed as part of this work displays the well-known phenomenon of overflow metabolism [42]. Reproducing this is possible using basic FBA with many arbitrary constraints, but is better captured using a model with properly tuned enzyme constraints, exemplified e.g. by the GECKO method [9]. The AutoPACMEN [14] package is streamlined to allow the addition and tuning of enzymatic constraints to a GEM. Its proof-of-concept, the iJO1366 � E. coli model, was therefore considered suitable for this work. Before tuning, all exchange reaction bounds were set to ±1000 mmol gCDW −1 h −1 . As specified in section GEM, COBRA and code availability, the Matlab program used for the simulation is detailed in S2 File, and the code can be found in S3 File.
In order to produce a dFBA scenario corresponding to the batch fermentation, the following tuning steps were performed: 1. As micronutrients are not supposed to be limiting in the given situation, exhaustion of the various ions and trace minerals was disabled (see S2 File for details on how, and S3 File FermentationSim/input/exchange_reactions/fermCL_ionsNotLimiting.csv for details on which).

2.
As stirring was coupled to DO in the fermentation, exhaustion of oxygen from the medium was also disabled.
3. The model-encoded uptake rates of glucose and oxygen were too high. These are closely linked, and as such both were lowered by limiting the oxygen uptake rate to a maximum of 15 mmol gCDW −1 h −1 .
4. The only secondary metabolite observed in the medium analysis of the batch fermentation to be excreted at an appreciable rate was acetate. However, the following were also excreted by the iJO1366 � model: lactate, dihydroxyacetone, ethanol, and pyruvate. Their excretion rates were therefore iteratively constrained to zero (both isomers where applicable), in the listed order.
5. Finally, the excretion rate of acetate was set to match observations. This resulted in a maximal acetate excretion of 2.2 mM gCDW −1 h −1 (see further discussion of this in the Results section).
Snapshots of the simulation for different parameter sets are showcased in Fig 2. The corresponding parameters are given in Table 1. The acetate excretion feedback is detailed in the Results and discussion section.
The medium recipe is given in S1 File. For the modeling, this was transcribed into a CSV file where the concentrations of the different constituent compounds are given in mmol L −1 . This can be found in S3 File. As implied by the tuning steps above, only the concentrations of glucose and ammonium were potentially limiting during the simulations. Of these, only glucose was exhausted.
After the above tuning steps, we applied the ecc constraint to observe the effect on the model's correspondence to experimental measurements. See Results and discussion.

Modeling a time-constrained change in enzyme composition
The central parameter in constraining the change in enzyme composition is implemented in the enzyme change constraint (ecc). This parameter limits the total change in enzyme composition with time, and is given in g gCDW −1 h −1 (referred to as "mass per hour", mph).
Before the ecc can be applied, a method for attaining the enzyme mass distribution is needed. In order to infer enzyme compositions, the flux through each reaction is multiplied by the corresponding bottom element of the sMOMENT stoichiometric matrix (S-matrix) and then multiplied by −1. The bottom element gives the molecular weight of the corresponding enzyme divided by the enzyme's turnover number. We get the total amount of mass currently dedicated to that enzyme per gCDW by multiplying this bottom element with the flux through the corresponding reaction. Multiplication by −1 corrects for the fact that the bottom row is all-negative in the sMOMENT formulation [14].

PLOS ONE
Any given enzyme-catalyzed flux v i is limited by the amount of enzyme g i dedicated to the enzyme, and its effective catalytic rate k cat,i . The amount of enzyme g i is in turn limited by the amount of mass that can be allocated to that enzyme, and the enzyme's molecular weight MW i . Expressed mathematically (similarly to Bekiaris & Klamt [14]), we get: where v i is the flux for reaction i in mmol gCDW −1 h −1 , MW i is the mass of the enzyme (complex) catalyzing reaction i in g mol −1 , k cat,i is the turnover number of enzyme (complex) i, and g i is enzyme "concentration" in mmol gCDW −1 .
If we assume the cell utilizes all its available enzyme in an optimal manner, we get: Using this equation, it is possible to keep track of the enzyme composition of the cell. Note that in this formulation, g i is purely implicit during computation. The mass for each enzyme is calculated directly from the reaction flux multiplied by the molecular weight of a given enzyme (complex) divided by its turnover number.
We apply an ecc rate limit by storing the enzyme composition as a vector. This vector has a length equal to the number of reactions, and specifies the amount of enzyme (in g gCDW −1 h −1 ), grams per gram of cell dry weight per hour) dedicated to each reaction. Subsequently, the ecFBA problem (without ecc constraints) is solved once, before the resulting enzyme composition is checked against the former composition. The difference in enzyme composition is calculated by multiplying the number of hours in the time step by the enzyme change limit, given in g gCDW −1 h −1 . If the difference is smaller than or equal to the amount of change permitted for that time step, the enzyme composition vector is updated to match the solution. If, however, the difference is greater than the limit for that time step, the enzyme composition is calculated by weighted linear interpolation. This interpolation occurs between the former enzyme composition and the optimal one, and a new ecFBA problem is solved with upper bounds on the flux values corresponding to this new interpolated enzyme composition.
Expressed mathematically, we have the enzyme composition vector of the current time step m t , giving the protein mass dedicated to each enzyme in gCDW: where E is the number of enzymes in the vector, which for the purposes of the sMOMENT model is inflated to the number of (uni-directional) fluxes in the model. For each time step, the composition vector of the former time step m t−1 is compared to the optimal composition vector m o acquired by solving a new ecFBA problem, thus giving the difference vector m d : If the sum of the positive elements in this difference vector is smaller than the enzyme composition change constraint γ for the current time step, then: Otherwise, the value for m t will be calculated by weighted linear combination. First we find the proportion of the change, α, permitted by the composition change constraint γ: Next, we acquire the new m t by weighted linear combination of m o and m t−1 : This option can be activated in the configuration file by setting the enzyme composition change limit for a given simulation. If set to −1, it allows the enzyme composition to change freely on every step without tracking, while set to any positive value (including 0), it constrains the change in enzyme composition to that number of gCDW h −1 . Any other number will limit the rate of change in enzyme composition to that fraction of total biomass per gCDW per hour. E.g. setting it 0.1 will allow 0.1 grams of enzyme to be replaced per gCDW per hour, which for the standard enzyme allocation setting of the sMOMENT model iJO1366 � of 0.0948 would amount to replacing 0:1 0:0948 ¼ 1:05, or 105%, of total enzyme per hour. To test the effect of constraining enzyme reallocation on a dFBA simulation, we first ran a model equivalent to standard dFBA, then added total enzyme constraints (decFBA), and finally added enzyme change constraints (decFBAecc). We selected a simple test scenario: a 1 L reactor with an initial biomass of 0.1 gCDW growing on 10 mM glucose. The bounds on glucose and oxygen uptake rates are specified in Table 2.
First, to get a regular dFBA equivalent, we used a version of iJO1366 � with the protein-pool constraint increased to allow arbitrarily high fluxes. This effectively inactivates the enzyme constraint, resulting in standard dFBA. Second, we ran a version of the same GEM where the protein-pool constraint remained at its default setting, thus performing standard decFBA. Third, and finally, we ran a version of the model where the rate of enzyme redistribution was constrained to a plausible rate value. This results in what we have called decFBAecc.

Experimental batch fermentation
During the exponential growth phases, growth rates were measured to 0.66 h −1 , while the terminal biomass values were measured to 8.2 gCDW L −1 . The experimental results, including sampling history and raw data from the fermentor's monitoring software and off-gas analyzer are given in the S4 File.
The fermentation data are plotted in Figs 2 and 3, where they are compared with different simulation results. As can be seen, the fermentation follows a typical timeline for diauxic growth of E. coli on glucose [7], with the depletion of glucose and concurrent accumulation of acetate, which subsequently is depleted. The concentrations of the various organic compounds in the medium were determined using the NMR protocol outlined in the Materials and methods section. Of note is the observation that biomass per volume starts to drop off after peaking, Table 2. Iterative change in constraints from dFBA to decFBAecc. The parameters that are subject to change are the total fraction of mass available for enzymes (ec; g gCDW −1 ) and the amount of enzyme in terms of total mass that can be replaced per hour (ecc; g gCDW −1 h −1 ).

PLOS ONE
likely due to starvation and cell death. Since the measurements of biomass are actual dry weight measurements, the drop-off is not a result of simple change in cell morphology. In order to get stronger signatures for the NMR medium analysis, high concentrations were used during the fermentation. As a side effect of this, the concentration of acetate seemed to reach saturation earlier than predicted by the computational simulations. This is potentially linked to thermodynamic feedback in the Pta-AckA pathway [43].

Applying enzyme change constraints to a dFBA model
In Fig 1 we show the effects of adding different constraints to the dFBA formulation, thus transitioning from standard dFBA to a decFBA approach, to a situation with enzyme composition change constraints (decFBAecc).
For these simulations, we set the O 2 maximum uptake rate to 15 mmol gCDW −1 h −1 and the glucose maximal uptake rate to 10.5 mmol gCDW −1 h −1 . The constraints are listed in Table 2.
The first parameter set ("dFBA", Table 2) corresponds to standard dFBA, realized through increasing the sMOMENT protein pool constraint from 0.0948 to 5000 g gCDW −1 of enzyme. Here, glucose is rapidly depleted as the amount of biomass grows and acetate is excreted, before the model fully depletes the glucose and then also rapidly depletes the acetate.
The second parameter set ("decFBA", Table 2) corresponds to decFBA as realized by running the sMOMENT approach unmodified in a dFBA setting. Here, glucose is consumed more slowly than for simple dFBA simulation, and acetate accumulates at a higher rate. This is typical when enzyme constraints are added, as incomplete metabolism of glucose through glycolysis produces more energy per time per amount of enzyme than does the full aerobic breakdown through TCA cycle and electron transport chain. Once glucose is depleted, the acetate concentration quickly follows suit at a similar rate to what is the result in the dFBA case.
The third parameter set ("decFBAecc", Table 2) corresponds to decFBAecc as realized by running the sMOMENT model unmodified in a dFBA setting while constraining the amount of enzyme mass that can be reallocated per time. The resulting curves seem identical to the decFBA case until glucose is depleted. This is expected since the only difference is the rate at which enzyme can be reallocated, and both decFBA and decFBAecc are initialized at optimal enzyme composition for the consumption of glucose. Following this, the acetate depletion rate builds up gradually, as the biomass curve under decFBAecc displays a lag phase. The final biomass concentration is lowered due to the application of these constraints reducing the total conversion efficiency of glucose to biomass.
Overall, we find that constraining the rate at which enzyme mass can be reallocated has a profound impact on the dynamics of cell metabolism. While the impact in terms of terminal biomass is small between decFBA and decFBAecc, the impact of the timing at which the terminal biomass is reached is significant. This has implications for microbial production of compounds where a change in growth conditions is often used as a trigger for production of a compound in question. In conclusion, the inertia in the cells' change in enzyme composition appears to warrant further investigation.

Modeling feedback in acetate excretion
When comparing fermentation data with model predictions, we found that the model overpredicted the external acetate concentration peak. Specifically, the acetate excretion seemed to slow with a higher external concentration of acetate. This phenomenon is described by Enjalbert et al. [43], and appears to be caused by thermodynamic feedback.
In order to capture this dynamic in our modeling framework, we tried to match their data to a number of functional forms using regression analysis: linear, polynomial (degree 2-5) ,   Fig 1. Transitioning from standard dFBA to decFBAecc. Shown in the figure are glucose (glc), acetate (ac) and biomass (DW) concentrations for dFBA, decFBA, and decFBAecc. Most reactions are unconstrained, but glucose uptake is limited to 10.5 mmol gCDW −1 h −1 , and oxygen uptake is limited to 15 mmol gCDW −1 h −1 . In the dFBA case, the iJO1366 � model has been adjusted with an arbitrarily large protein pool constraint (5000 g gCDW −1 of enzyme, which allows every flux to reach 1000 mmol gCDW −1 h −1 ) to produce standard dFBA results; glucose is rapidly depleted, while acetate is accumulated before being rapidly depleted as well. In the decFBA case, the iJO1366 � model is simulated without modifying the protein pool constraint (0.0948 g gCDW −1 of enzyme), allowing slightly slower growth; glucose is depleted slower than in the dFBA case, and more acetate accumulates before being depleted at a similar rate. In the decFBAecc case, the protein pool constraint is the same as default setting, but the rate at which the enzyme composition is allowed to change has been limited to 1% per hour (as opposed to being instantaneous); the trajectories are identical until the metabolism shifts from growth on glucose to growth on acetate, at which point there is a short lag phase. https://doi.org/10.1371/journal.pone.0280077.g001

PLOS ONE
and logarithmic. We found that a logarithmic function provided an excellent fit, but that the parameters from our regression model did not allow our model to reach the experimentally measured concentration of external acetate.
The difference could owe to differences in the strains of E. coli used, or due to differences in medium or exact growth environment. The expression we devised to fit our data is given in Eq 8: where Ex Ac,max is used as the upper bound on acetate excretion, and c Ac is the concentration of acetate in the environment. V Ac,max = 2 is the max excretion rate for acetate, while c Ac,s = 13 is the external acetate concentration at saturation.  Table 1.

Fitting model parameters to batch fermentation
The base dFBA computational modeling, as seen in Fig 2 panel A, uses the iJO1366 model. It depletes glucose too rapidly and produces an acetate peak that is too high. Terminal biomass is not too far off, however. In Fig 2 panel B, oxygen uptake has been limited to a more realistic value (see Table 1), which caused excretion of lactate and dihydroxyacetate. As these secondary metabolites were not observed to be excreted in our batch fermentation, they have been blocked. Overall, this results in a slower glucose depletion and a higher acetate peak, without an appreciable change in terminal biomass. In Fig 2 panel C, acetate excretion and uptake are limited to more realistic values (see Table 1) [43]. This caused excretion of pyruvate, which was subsequently also blocked, as no appreciable amounts were detected in the medium. This step appears to make growth, glucose, and acetate curves align much better with experimentally measured values. Finally, in Fig 2 panel D, we added an expression for negative feedback from external acetate concentration on acetate export (see Modeling feedback in acetate excretion for details) motivated by the results by Enjalbert et al. [43]. The resulting concentration curve for external acetate appears to align closely with our experimentally measured values, except its utilization is too rapid after glucose depletion.
After all tuning steps are done, the growth rate and carbon source consumption on glucose are 0.65 h −1 and 7.9 mmol gCDW −1 h −1 respectively, while on acetate (after the diauxic shift) they are 0.25 h −1 and 10 mmol gCDW −1 h −1 .

Consequences of the enzyme change constraint
Finally, we investigated the effect on the diauxic shift, from glucose to acetate, of applying a constraint on the rate of change in the cell's enzyme composition, as formulated in the decF-BAecc method (see Materials and methods). Here, one would expect to see a slower adaptation than for dFBA/decFBA without the ecc constraint, as time is required to reallocate part of the enzyme pool to allow optimal utilization of acetate rather than glucose.
The resulting simulation outputs are plotted in Fig 3 for a range of sample values for the ecc parameter. The consequences of applying the constraint on the rate of enzyme reallocation are most apparent in the depletion of external acetate, and, to a lesser extent, in the terminal biomass.
The more constrained the reallocation of enzymatic mass becomes, the slower the rate of acetate depletion approaches its optimum. This has a consequence of lowering the terminal biomass, as more acetate is consumed for energy by the non-growth associated maintenance sink function before it can be used as an energy and carbon source to fuel growth. This is to be expected, as the optimal enzyme distribution for utilization of glucose is not the same as the optimal enzyme distribution for utilization of acetate. As can be seen, several export reactions need to be blocked, and several rates limited in order to reproduce realistic behavior, even with the total enzyme pool constraints. The order and nature of adjustments are listed in Table

PLOS ONE
In this particular case, a rate of replacing 2% of total biomass worth of enzymes per hour produces the best fit with experimental data according to our analyses, with a root mean square deviation (RMSD) of 1.04. The entire interval 1.8-6.1% for the value of the ecc appeared to fit experimental data quite well, with an RMSD below 1.10. The rate is given as percentage of total cell mass (not enzyme mass) per hour (denoted mph in the legend) that can be reallocated. The rate of change of enzyme composition has a dramatic effect on the transient behavior of the model during diauxic growth, and a rate of 2% mph appears to best match the data measured from the fermentation experiment. The difference in terminal biomass is caused by the non-growth-associated maintenance (NGAM) function, which forces hydrolysis of ATP at a rate of 3.15 mmol gCDW −1 h −1 . https://doi.org/10.1371/journal.pone.0280077.g003 It should be noted that the ecc rate is likely to be highly dependent on rate-limiting components in the enzyme reallocation process, i.e. gene expression, which should place a lower bound on how rapidly this can occur. A more accurate model would also have the ecc rate represented as a vector of separate values for different reactions, though a lack of appropriate data makes this infeasible at present.
With the general approach used for decFBAecc, one would expect a similar effect for diauxic shifts involving other substrates, such as glucose/glycerol or glucose/lactose. In their paper from 2006, Bettenbrock et al. modeled diauxic shifts on a number of substrate combinations, with a focus on the signal transduction processes involved [44]. This model, in their own terms, "describes the expression 17 key enzymes, 38 enzymatic reactions, and the dynamic behavior of more than 50 metabolites" [44]. While their model by and large yields good fits with the data, it is custom-built to describe the process, and may leave out other elements of interest when modeling the process. The strength of the decFBAecc approach lies in its foundation, the robust framework and formalisms surrounding GEMs, FBA and constraint-based analysis. This foundation allows standardized approaches to applying enzymatic data, such as in the case of sMOMENT [14]. This in turn allows the representation of enzyme reallocation dynamics with minimal additional overhead and parameters using decFBAecc.

Conclusion
In this work, we have presented a computationally robust and flexible framework for simulating fermentation with dynamic flux balance analysis. We use this framework to make iterative changes in computational constraints applied to a simple example of diauxic growth on glucose. This approach illustrates the consequences of adding constraints on the rate of reallocation of enzyme. Based on these findings, we perform a simple fermentation experiment for E. coli as a case study, which we attempt to recreate using our computational simulation framework.
We have two primary findings. First, we find that applying constraints on the rate of enzyme mass reallocation, in addition to applying constraints on total mass allocation to enzymes, produces markedly different modeling results for key processes of interest. The new constraints have a sound rationale, and allow the model in this case to better capture experimental observations. One key consequence is that they markedly change the timing at which maximum biomass is reached during diauxic growth. They do not, however, markedly change the terminal amount of biomass in the considered case, meaning their utility in verifying a BOF may be limited. Another key consequence is that they drastically alter the rates and timings of substrate depletion during a change in growth environment, such as a diauxic shift. This could have important implications for industrial production of added-value compounds. And as the method is relatively straightforward to apply, at least for organisms for which highquality GEMs exist, it can be readily employed for such purposes.
Second, we strengthen the notion that non-linear feedback dynamics on acetate excretion are necessary to capture the behavior of E. coli in environments with high initial concentrations of glucose. This agrees with the findings of Enjalbert et al. [43], and has important consequences for modeling of metabolism at high external concentrations of export products.  Almaas.